Fast calculation of the electrostatic potential in ionic crystals by direct summation 

method 
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An efficient real space method is derived for the evaluation of the Madelung's potential of ionic 
crystals. The proposed method is an extension of the Evjen's method. It takes advantage of a general 
analysis for the potential convergence in real space. Indeed, we show that the series convergence is 
exponential as a function of the number of annulled multipolar moments in the unit cell. The method 
proposed in this work reaches such an exponential convergence rate. Its efficiency is comparable to 
the Ewald's method, however unlike the latter, it uses only simple algebraic functions. 

PACS numbers: 



I. INTRODUCTION 

Since 90 years a number of methods have been pro- 
posed to calculate the electrostatic potential in ionic crys- 
tals. These methods can be separated into two categories, 
the direct summation methods and the indirect summa- 
tion ones. The former uses a real space summation of 
the electrostatic potential generated by the ions within a 
finite volume {£). However, when enlarging the volume 
£, such partial summations are conditionally convergent. 
The convergence depends on the specific shape oi £. In 
addition, when achieved, the convergence is quite slow. 
The indirect summation methods do not present these 
drawbacks since the long range part of the potential is 
calculated in the reciprocal space. Indeed, the summa- 
tion is divided into two parts, a short range one, eval- 
uated by a direct summation in real space and a long 
range one evaluated in the reciprocal space. Among 
these methods, the most widely used is the Ewald's 
mcthodi, which is actually considered as the reference 
for Madelung potential calculations. 

Despite its quality the Ewald's method is not easily 
usable in different domains of physics. This is for in- 
stance the case in clusters ab-initio calculations used for 
the treatment of strongly correlated systems, the study of 
diluted defects in materials or adsorbates or for QM/MM 
type of calculations. For this types of calculations, real 
space direct summation methods are used. There is thus 
a need for efficient and accurate techniques for the deter- 
mination of the Madelung potential in real space. 

The convergence problems found in real space summa- 
tion are linked to the shape of the summation volume, 
£, and more specifically the charges at its surface. In 
order to insure the convergence of the summation, the 
surface charges are renormalized. Several methods have 
been proposed for this purpose. 

The most common and simple one is the Evjen's 
method^. This method uses a volume £, built from a 
finite number of crystal unit cells, and renormalizes the 
surface charges by a factor 1/2, 1/4 or 1/8 according 
whether the charge belong to a face, edge or corner of £. 
This method insures, in most cases, the convergence of 



the electrostatic potential when £ increases. However, in 
some cases such as the famous CsCl it does not converge 
to the proper valuei^ 

Other authors^ii proposed to renormalize not only the 
surface charges, but also the charges included in a thin 
skin volume. The adjustment of the renormalization fac- 
tors are, in this case, numerically determined so that to 
reproduce the exact potential, previously computed us- 
ing the Ewald's method at a chosen set of positions. Such 
a method presents the advantage of reaching a very good 
precision. However several drawbacks can be pointed 
out : i) the previous calculation of the electrostatic po- 
tential using the Ewald's method at a large number of 
space positions, ii) the necessity to invert a large linear 
system to determine the renormalization factors and iii) 
finally the fact that the latter are not chosen on physi- 
cal criteria. Indeed, this last point induces the possibility 
that the renormalization factors can be either larger than 
one or negative. It results that even if the electrostatic 
potential is very accurate at the chosen reference posi- 
tions, its spatial variations can be unphysical and thus, 
the precision can strongly vary when leaving the reference 
points. 

Marathe et afe suggested a physical criterion, based on 
the analysis of the convergence of the real space summa- 
tion, for the choice of the renormalization factors. In- 
deed, it is known that the direct space summation con- 
verges to the proper limit if the volume £ presents null 
dipole and quadrupole momenta^. The authors of ref- 
erence |5| showed, on the simple example of a linear al- 
ternated chain, that it is possible to find a finite number 
of charge renormalization factors allowing the cancella- 
tion of these two multipolar moments. They also assert 
that the cancellation of additional multipolar moments 
increases the speed of convergence. Unfortunately they 
did not prove this affirmation and more importantly, they 
did not proposed a practical way to determine the renor- 
malization factors in order to reach this goal. 

In the present paper we propose a systematic method 
for the determination of the renormalization factors al- 
lowing the cancellation of a given number of multipo- 
lar moments as well as a careful analysis of the direct 



space summation convergence as a function of the num- 
ber of canceled multipolar moments. The next section 
will present the convergence proof, section 3 will develop 
the method for determination of the renormalization fac- 
tors and section 4 will present the optimization of the 
method and illustration on a typical example. 



II. CONVERGENCE ANALYSIS 
A. Potential at a point 

As already mentioned in the introduction, several pa- 
pers already exist on this subject. However the results 
are only partial and there is not complete analysis of the 
convergence issue. We will thus present in this section a 
global analysis of the electrostatic potential convergence 
in a real space approach and an estimation of the error. 

We want to evaluate the limit of the following scries 
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where J? is a vector of the Bravais's lattice, j refers to 
a charge qj located at the position tj of the unit cell C. 
{Sm n G N} is a set of volumes such that 
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For the sake of simplicity we require that the set of £„ 
also presents the following conditions 
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The well known problem of this series is that the limit 
depends on the particular choice of the unit cell C and of 
the volumes £„■ In fact different shapes of the charge set 
will result in different limits due to the surface effects. 
However, it has been demonstrated that this conditional 
convergence disappears if one considers a unit cell with 
zero dipolar and quadrupolar moments^. We will con- 
sider in the following that the cell C fulfills this condi- 
tion. In this case, one only obtains the so-called "bulk 
contribution" as in the Ewald's and related methods. 

The error AVn{fo) on the electrostatic potential eval- 
uation, Vn{r()), can be written as 
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For large values of n, AVn{fo) can be evaluated by a 
multipolar expansion. For practical reasons, we will use 
an expansion expressed in spherical coordinates. Indeed, 
for a given order, this expansion contains less terms than 
the usual multipolar expansion based on Cartesian co- 
ordinates. The use of the later multipolar expansion is 



still possible but is more complex (see ref@). In spherical 
coordinates, the multipolar expansion of the error made 
on Vn (ro) reads^ : 
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where Aiimi'ro) are the multipolar moments of unit cell 
C at fb- They can be expressed as 
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(i?, 6, (p) are the spherical coordinates of the Bravais's 
vector R and {pj,9j,(j)j) the spherical coordinates of 
rj —fo. The ^~™ are Schmidt semi- normalized spherical 
harmonics : 



Yr"'{o„<i>,) = y fq:^^r(cos0)e™^ (s) 

where P™ are Legendre functions. 

In order to overvalue the error, one needs to overvalue 
the spherical harmonics. We thus consider the addition 
formula : 
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where Pi is a Legendre Polynomial and 7 is the angle 
between {R, 6, (p) and (i?', 9' , (f>'). It comes for 9^9' and 
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Using this result, one obtains the following overvalue 
for the moments : 



\Mim\ <Q(ro + a)' 



(12) 



where a is the typical size of C (i.e. the diameter of 
the circumsphcrc of C) and Q is the sum of the absolute 
values of its charges 
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One should notice at this point that Yi^{9, (p) has the 
same parity as I. The contributions of R and —R unit 
cells to AVn {ro) thus cancel when / is odd. Using equa- 
tions [11] and [121 one gets the following overvaluation : 
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where 2p is the order of the first even, non-zero moment. 
Let us now overvalue the sum over the 1/R powers 
by a volume integral. For each cell C{R) located at R, 
the norm of the position vectors r belonging to the vol- 
ume C{R) is smaller than R + a. One can thus overvalue 
l/i?2fc+i by 
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LJ being the volume of the unit cell C. If 7^„ is the radius 
of the insphere of £„, it comes 
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The later sum converges for TZn large enough, i.e. for 
T^n > ''o + 2a. The decreasing function (4fc + l)/(2fc — 
2) can be overvalued by its value in k = p. Further 
summation over k leads to the following expression 
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with 
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The electrostatic potential at rg thus converges as l/7^Jj ^ 
where I is the first, even, non-zero moment of the unit 
cell C. 



B. Difference of potential betw^een two points 



In several applications, as for instance in cluster ab 
initio calculation, the problem depends on the spatial 
variations of the potential and not on its absolute value. 
In such cases it is sufficient to cancel the dipolar moment 
of the unit cell in order to ensure the convergence of the 
calculation. The convergence rate can also be expected 
to be faster than for the calculation of the potential at a 
point as we will show in this section. 

Let us overvalue the error made on the calculation of 
a difference of potential between two points located at 
fo + fi and fo — fi : 
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where {R-,0 -,(/>-) and (_/?+, 6'+, 0+) are the spherical 
coordinates oi R — fi and R + fi respectively. As in the 
preceding section (i?, 9, (j)) will be the spherical coordi- 
nates of R and (ri ,9i,(j>i) those of fi . 

In order to express the previous expression as a func- 
tion of 1/R, we use following expansion of solid spheri- 
cal harmonics (for simple derivation see ref. [lo, see also 
ref . [nllir 
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where the sum over m spans all integer values. Neverthe- 



less, only a finite number of terms will contribute, since 
P™ = if |77i| > L Setting I = L + n and introducing 
spherical harmonics leads to : 
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Considering P+ in this equation instead of P_ is equiv- 
alent to the transformation 
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that results in an overall (—1)'"^ factor. Inserting rela- 
tion [21] into eq. HD] and inverting the summation over / 
and L leads to the expansion : 
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Considering the parity of Y/,„ , one can see from eq. [23] 
that the contributions from cells located at R and ~R 
cancel when I is odd. Moreover, due to the last term of 
eq. [211 the Aim coefficients are zero when I ~ L is even. 
Only terms with I even and L odd have a non zero con- 
tribution, thus only moments MLAiiT^o) with odd order 
will contribute to the error. The consequence is that the 
first non-zero contribution in equation 1231 corresponds to 



I = 2p+2 where 2p-|-l is the first, non-zero, odd moment 
of the unit cell. 

Let us now find an overvalue of the Aim terms. It is 
easy to show using a recurrence relation on the values of 
TO and M, that if |77i| < I, \M\ < L and \M-m\<l-L, 
the following relation holds : 
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Using the previous overvaluation of the moments (cq. [T 
one obtains : 
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where the summation runs only up to / — 1 since / and L 
are of different parity. It comes 

\Aim\ < 2Qi2l~l)l{ro + a + ny-^n (27) 

As in previous section, the sum over TZ can be overval- 
ued by a volume integral (cf . eq I16p . Overvaluation of 
the error thus reads : 
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The sum over k converges if 7^„ is larger than ro + 2a + 
ri. It can be calculated using derivative of power series. 
After simplification, one obtains : 
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As one increase the size of the set of charges, the dif- 
ference of electrostatic potential between two points con- 
verges like l/7^„ ^, where I is now the first, odd, non-zero 
moment of the unit cell C. This convergence is slightly 
faster than the convergence of the absolute value of po- 
tential. When the order of the first non zero moment 
is even, the convergence rates differ by a factor 1/R^, 
otherwise they are similar. 



C. Electric field at a point 

The convergence problem of the electric field at a point 
is very similar to the problem of the difference of poten- 



tial between two points. However since it could be of 
practical interest, for instance for molecular dynamists 
in the calculation of ionic forces, we will provide in this 
section the analysis of the electric field convergence. 

The error on the evaluation of the electric field at a 
given point rg is related to the error of the potential 
difference between two nearby points as 
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where E" is the a component of the electric field and Uq 
is the unit vector in the a direction. 

AE"{ro) can thus be overvalued using equation]^ 
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As expected, one sees that the electric field converges 
with the same rate as the potential energy difference be- 
tween two points, that is as l/7?.Jj~^, where I is the first, 
odd, non-zero moment of the unit cell C. 



III. PARTIAL CHARGES 

As depicted in the previous section, convergence can 
be considerably increased if one cancels several multipo- 
lar moments of the unit cell. In general the Evjen method 
allows to only cancel the dipolar moment, and thus pro- 
vides a convergence of the potential differences in 1/7?.^. 
In order to really take advantage of the former property, 
one needs a method allowing the cancellation of several 
multipolar moments. 

In this section we will establish a method to con- 
struct unit cells with a chosen number of zero multipolar 
moments. The method, based on the usage of partial 
charges, is general and can be applied to any Bravais's 
crystal. 

Let {So, bo, Co) be the lattice vectors of the Bravais's 
crystal, and Co the associated unit cell. In order to in- 
troduce partial charges, we consider a larger cell Ci of 
dimensions (l x ao,l x bo J x co), that we will refer as 
the "construction ceW . The construction cell thus con- 
tains l'^ original unit cells Co which positions in C; can be 
labeled by p, q and r indices, ranging from 1 to I. 

If we note Uc the number of charges qi in the orig- 
inal cell Co, the cell C; now contains n^ x P charges. 
These charges will be corrected by a factor XL (where 
i refers to the charge qi). When on rebuild the lattice 
using the construction cells Cj, the cells overlap, and the 
final charge at position fi corresponds to the superposi- 
tion of partial charges from several construction cells. It 
is straightforward to show that the condition to retrieve 
the nominal value of the charges qi reads : 
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At this stage, considering the latter Uc equations, the 
cell Ci contains nc{l^ — 1) free parameters that could be 
used to cancel multipolar moments. For the sake of sim- 
plicity and generality (i.e. for the method not to depend 
on the particularity of a given crystal), we will impose 
further conditions on the A'^ coefficients. 

We first reduce the problem to a one dimensional prob- 
lem by setting : 
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where the three coefficients A"'' 



A^'' and A^'' are used to 



cancel multipolar moments of the ID problems obtained 
when the cell Co is respectively projected on the three 
axes of the crystal. For each one dimensional problem, 
the construction cell contains I projected unit cells. The 
condition on the coefficients now reads : 
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It is easy to show that, if these coefficients cancel a fixed 
number of multipolar moments in each one dimensional 



problems, the A*^ coefficients will cancel the moments 
of same order in the original three dimensional problem. 
We further impose to the coefficients to only depend 
on the fractional coordinates (cti, Pi,ji) of the charges qi, 
in the corresponding direction: 



Xp = Xp{ai) Xq{(3i) Ar(7i) 



(36) 



The Xp{x) functions arc thus the same for the three di- 
rections and for all charges. As a consequence their ex- 
pression is the same for all crystals. For a given value of 
X, and considering the condition for the reconstruction of 
the crystal (eq. [55]) , wc are left with I — 1 degrees of free- 
dom. We thus impose to the Xp{x) functions to cancel 
/ — 1 multipolar moments. This can be done by setting 
the moments created at the center of the construction 
cell by a unique charge q : 



9"o ^K{x) yx+p-\ 
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= quo mi^k (37) 
(0 < fc < / - 1) 



where uq = ao,bo, or co, x is the fractional coordinate 
of the charge and mi^k are constant values (independant 
of the crystal specifications). The equation obtained for 
k = corresponds to the condition for the reconstruction 
of the crystal (m;^o = !)■ The moments of the construc- 
tion cell, can thus be obtained by summing the contri- 
butions of all charges. As the unit cell is neutral, these 
contributions cancel out. 

The equations [351 and [37] thus define sets of partial 
charges that allow to construct cells with I — 1 zero mul- 
tipolar moments. The shape of these partial charges de- 
pends on the choice of the mii; constants values. In order 
to find the most reasonable choice of partial charges, we 
search for mi_k constants that satisfy the following phys- 
ical conditions : 

1. Vp, Xp{x)e[0,l]. 

2. The partial charges vary continuously, i.e. 
\fp, Xp{x) arc continuous and Ap(l) = Ap+i(0). 

3. Vp, Xp{x) decreases monotonously when moving 
away from the center of the construction cell. 

4. The values of the mi^k are as small as possible. 

The latter condition ensures that the partial charges are 
larger in the center of the construction cell and smaller on 
its edges, and hence that this cell is close to the original 
one. 

Wc did not find a way to derive the solution of this 
problem in a general way, for any value of I. We thus 
determined the solution for fixed values of / up to / = 6. 
In all these cases the first Ap function presents the same 
shape : 
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We reasonably assume that this expression is vahd for 
any value of /. As we will see, it is possible to show, 
a posteriori, that the Xp functions fulfill the first three 
conditions. 

We will now determine the function Aj, for any I, using 
the above expression of Ai and relation 1371 In equa- 
tion [371 the TTifc.i constants can be replaced by the value 
of the moments obtained for x = : 

p=i ^ ^ p=i ^ 

(39) 
with < k < I — 1. The matrix of this linear system is 
the transpose of a Vandermonde Matrix. Inversion of the 
system leads to : 



These coefficients are represented on fig. [2l The abscise 
values have been shifted by {I — l)/2, so that the origin 
corresponds to the position of the effective edge of the 
fragment (i.e. the position of the edge obtained when us- 
ing n original cells Cq without partial charges). One can 
see that the renormalization of the charges is relatively 
small. Indeed, even in the case I = 10, this renormal- 
ization is weaker then 1% for charges at distances larger 
than 2uo, (wq — ^o , bo or cq) . 

As already mentioned, the Ap and fip are polynomial 
functions of order l — l. It is also interesting to notice that 
at the junction point between the ^ip{x) (resp. Xp{x)) and 
/ip+i(a;) (resp. Ap+i(a;)) functions, the renormalization 
function and all its the derivatives are continuous, except 
for the last one {[I — 1]*'' derivative). 
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The Ap(0) coefficients can now be obtained from the 
expression of Ai. Let us consider the latter equation eval- 
uated for p = 1, and for / integer values of x : 
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Replacing Ai(n) by its expression and inverting this lin- 
ear system leads to the expression of Ag(0) coefficients : 
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Finally, using eq. [40l and l42l one obtains the expression 
of the Xp{x) which are polynomial functions of order / — 1 
in X (see fig. [1]) : 
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.(43) 
These functions are segments of the uniform sum distri- 
bution , i.e. the distribution Pi{x) of the sum of/ uniform 
variates on the interval [0, 1] : 

Xp{x) = Pi{x+p- 1) (0 < X < 1 and 1 <p< /- 1) 

(44) 
From this relation, it is obvious that the Xp (x) functions 
satisfy the first three conditions mentioned above. 

We now consider a fragment of crystal made of n con- 
struction cells Ci . As these cells are composed of I original 
cells Co, they partially overlap, and the size of the frag- 
ment corresponds to n + / — 1 cells Cq. n' = n — I + 1 
cells Cq in the center contains charges with the nominal 
values qi , and I — 1 cells Cq on each side of the fragment 
contains partial charges. The latter partial charges are 
proportional to the coefficients : 

p 
//p(x)==^A,(x) (l<p</-l) (45) 

9=1 



IV. OPTIMIZED METHOD 

We first use the cells Ci to illustrate the convergence 
of the standard, real-space calculation of the potential 
established in section |TT1 In order to observe the general 
behavior of the different methods, we chose the a quartz 
structure, that possesses a reasonable number of atoms 
per unit cell and not too many symmetries. A cell Ci is 
constructed for a fixed number of p. The cells are used to 
produce sets of charges of increasing size. The potential 
and the electric field are calculated at the position of the 
twelve atoms of the central unit cell Cq. 

Fig.[3]a) represents the maximum of the error made on 
these potential values as a function of the width of the 
set of charges, fig. |3] b) represents the maximum of the 
error made on the sixty six differences of potential, and 
finally fig.[3]c) represents the maximum of the error made 
on the electric field. As expected the errors decrease as 
power functions of the width D of the set of charges. It 
fully agrees with eq. [TSl [2H] and [321 In particular, the 
fact that the cancellation of a moment of odd order do 
not increase the convergence rate of the potential at a 
point, clearly appears. Similarly, the cancellation of even 
order moment do not improve the convergence rate of the 
calculation of differences of potential and electric field. 

One can see from the previous figures that this stan- 
dard approach is not the more efficient. Indeed the in- 
crease of the number of zero multipolar moments clearly 
yield a faster convergence rate than the increase of the 
volume of the system for a fixed value of I. Let us there- 
fore fix the width n' of the volume containing the nominal 
charge, and let / increase. The variation of the maximum 
error made on the potential, on the potential differences 
and on the electric field are respectively represented on 
fig- HI a), b) and c). One sees that the present approach 
leads to an exponential convergence of the potential in all 
cases. The convergence is very fast, since an increase of 
the set of charges width by two unit cells results in a pre- 
cision increase by a factor better than ~ 10. Increasing 
the number n' of central cells without partial charges has 
a small influence on the convergence speed. For n' > 3 
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FIG. 1: (color online) Ap(a;) functions obtained for values of I ranging between 2 and 10. Intervals correspond to cells Co that 
compose the construction cell Cj. In each intervals, fractional coordinates x range between and 1. 
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FIG. 2: (color online) Renormalization of the charges at the "left" edge of a crystal fragment, obtained for values of / ranging 
between 2 and 10. The origin of the abscises corresponds to the position of the edge obtained when using the original cell Co 
(dotted line). 



the method even becomes less efficient since increasing 
n' increases the size of the total set of charges. The best 
convergence is obtained for n' = 3. It corresponds to the 
case where the cell in which the potential is calculated is 
surrounded by one shell of cells with the nominal charge 
values. 

Finally we compare our method to the famous Ewald's 
method which mixes calculation in real space and recip- 
rocal space. This method introduces Gaussian distribu- 
tion of charge ea;p(— a^r^), where the a coefficient can 
be adjusted. Increasing a coefficient increases the con- 
vergence rate of the real space sum, but slows down the 
sum in reciprocal space. A width of Gaussian propor- 
tional to the characteristic length of the unit cell, which 



corresponds to a^ 



'1/3 



is generally assumed to give 



a good compromise. We calculated the error made on the 
value of the potential using the Ewald's method for dif- 
ferent values of a around ao ■ The results arc represented 
on figure [51 as well as the error of our method obtained 
for n' ~ 3. Figure [5] reports the error on the potential 
as a function of the number of construction cells used in 
the calculation. Let us point out that, while this variable 
is pertinent for the global convergence rate analysis, for 
each charge, the Ewald's method requires an error func- 
tion evaluation resulting in a non negligible pre- factor, 
not present in our method and not taken into account 
in figure [5] One sees that the convergence rate of the 



present method is comparable with the Ewald's method. 
If one is only interested in the potential evaluation at 
a single point, the Ewald's method with an optimal a 
parameter is somewhat faster than the present one. One 
the other hand, once the renormalization have been com- 
puted, the value of the potential at any other point of the 
of the central area can be calculated with a similar pre- 
cision at little cost. More important, properties using 
potential integrals or complex potential functions can be 
more easily evaluated since our method used only alge- 
braic functions. 



V. CONCLUSION 

Number of authors have searched for a fast converging 
method for the evaluation of the electrostatic potential 
in real space. Similarly, many works where done yield- 
ing partial results on the convergence rate of such real 
series. The present work fills the gaps and proposes a 
general analysis of both the convergence of the poten- 
tial at one point and of the convergence of differences of 
potential. Indeed, we gave a general and rigorous proof 
of the relation (claimed by other authors) between the 
power law convergence of the series and the number of 
zero multipolar moments of the crystal construction cell. 

Based on these convergence analyses we derived a gen- 
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FIG. 3: (color online) a) Maximum of error made on the 
calculation of the potential at the position of the atoms of 
the central unit cell. D = n' -\- 21 — 2 is the width of the 
sets of charges in each crystallographic direction, where n' 
is the width of the central zone where the charges have their 
nominal values. / is the order of the first non zero moments of 
the cell Ci used to construct the set of charges, b) Maximum 
of error made on the calculation of the difference of potential 
and c) on the electric field. 



FIG. 4: (color online) a) Maximum of error made on the 
calculation of the potential at the position of the atoms of 
the central unit cell. The abscise corresponds to the width 
D = n' + 21 — 2 oi the sets of charge in each crystallographic 
direction, where n' is the width of the central zone with nom- 
inal charge values, b) Maximum of error made on the cal- 
culation of the difference of potential and c) on the electric 
field. 



eral real space method with an exponential convergence 
rate, comparable with the Ewald's method. The expo- 
nential convergence is reached as a function of the num- 
ber of canceled multipolar moments in the construction 
cell. The crystal is indeed constructed using overlapping 
construction cells with rcnormalizcd charges. Wc derived 



a general analytical expression of the renormalization fac- 
tors, for any given number of zero multipolar moments. 
Finally, we would like to point out that our method 
warrants continuous and smooth variations of the renor- 
malization factors. This property is of particular interest 
for molecular dynamic usage since it insures continuous 
and smooth variations of the ionic forces as a particle 
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crosses the cell boundaries. One can see the present /i 
functions as optimized cut-ofF functions. 



FIG. 5: (color online) Error on the Madelung's potential 
evaluation using the present method (bold black solid curve) 
and the Ewald's method with different parameters a (colored 
dashed curves). 
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